module chkIS
	use analyticalIS
	use compute
	use maxRPmod
	use ssmaster
	use CImod
	implicit none
		
	contains
	
	subroutine chkCIs
		integer	::	lb1, lb2,l,i,cc
		real	::	Y1(k),Y2(k),Zsum
		integer, parameter	:: ngrid=300
		integer	::	mugrid(ngrid),mug2(ngrid)
		real	::	mu, muhat,s,mu2,mumax,mumin, ci(2)
		real, allocatable	:: bys(:,:)
		
		allocate(bys(2*k+1,0))
		cc=0
		do lb1=1,nblocks
			lb2=lb1+1
			if(lb2>nblocks) lb2=1
			do l=1,npb
				Y1=Y0s(1:k,l,lb1)
				Y2=Y0s(1:k,l,lb2)
				if(Y1(1)==maxconsttail .or. Y2(1)==maxconsttail) cycle
				cc=cc+1
				Zsum=Y0s(0,l,lb1)-Y0s(0,l,lb2)
				if(Zsum>=(n_sample-2*k)*Y1(k) .or. -Zsum>=(n_sample-2*k)*Y2(k) ) cycle

				ci=[getCIlb(y1,y2,Zsum,sstst),-getCIlb(y2,y1,-Zsum,sstst)]
				
				muhat=(sum(Y1)-sum(Y2)+Zsum)/n_sample
				s=5*sqrt(1+sum((Y1-muhat)**2)+sum((Y2+muhat)**2))/real(n_sample)
				mu=muhat
				do 
					mu=mu+s
					if(reject([Zsum-(n_sample-2*k)*mu,Y1-mu,Y2+mu],sstst)) exit
				enddo
				mumax=mu
				mu=muhat
				do 
					mu=mu-s
					if(reject([Zsum-(n_sample-2*k)*mu,Y1-mu,Y2+mu],sstst)) exit
				enddo
				mumin=mu
				do i=1,ngrid
					mu=mumin+(mumax-mumin)*(i-1.0)/(ngrid-1)
					mugrid(i)=boole(reject([Zsum-(n_sample-2*k)*mu,Y1-mu,Y2+mu],sstst))
					mug2(i)=boole(mu<ci(1) .or. mu>ci(2))
				enddo
				if(.not. all(mugrid==mug2)) then
					call mydispivec(mugrid)
					call mydispivec(mug2)
					print *,""
				endif
				if(sum(abs(mugrid(2:)-mugrid(1:ngrid-1)))==2 .and. sum(abs(mugrid-mug2))>2) then
					call mydispivec(mugrid)
					call mdisp([Zsum,Y1,Y2])
					bys=reshape([bys,[Zsum,Y1,Y2]],[2*k+1,size(bys,2)+1])
					call savemat(trim(folder)//"bys.txt",bys)
					call mdisp(real([size(bys,2),cc]))
				endif
			enddo
		enddo
		call printtime
		print *,"done with CI chk"
		stop
		
	end subroutine
end module
	
program mfort
!dir$ nooptimize
	use analyticalIS
	use compute
	use maxRPmod
	use ssmaster
	use setlammod
	use derivs
	use chkIS
	USE omp_lib
    implicit none
	
	integer		:: i,im,j
	integer	:: nlist(3)=[50,100,500],in
	
	real	:: levellist(23)=		[0.05,0.01,0.1,		0.002,	0.004,	0.006,	0.008,		0.02,	0.03,	0.04,		0.06,	0.07,	0.08,	0.09,	0.12,	0.14,	0.16,	0.18,	0.2,	0.25,	0.30,	0.40,	0.50]
	real	:: levelfudgelist(23)=	[0.002,0.001,0.005,	0.0002,0.0002,	0.0005,	0.0005,		0.002,	0.002,	0.002,		0.003,	0.003,	0.003,	0.005,	0.005,	0.005,	0.005,	0.005,	0.005,	0.01,	0.01,	0.01,	0.01]

! comment in if only test of level 5% and 1% are to be determined, such as for k=12
	!real	:: levellist(2)=		[0.05,0.01]
	!real	:: levelfudgelist(2)=	[0.002,0.001]

	integer	:: il
	
	call inittime
	
	call rnset(25)
!$omp parallel do
	do i=1,50
		call rnset(22+134753*omp_get_thread_num())
	enddo
	
	call mkGQxw(GQxw)
	call setmupper
	call savemat(trim(folder)//"GQxw.txt",GQxw)

	print *,"k=",k
	print *,"n0",n0
	print *,"varestfac=",varestfac
	print *,"swsmall=",swsmall, "swflat=",swflat
	print *,"xsi null range"
	call mdisp([xsimin,xsimax])
	print *,"xsi alternative range"
	call mdisp([xsimin1,xsimax1])
	
	stage=1

	call setanaIS
 	call prep

	do il=1,size(levellist)	

	if(tstconst) then
		if(il==2) then
			stage=2
			level=0.05
			call loadlam
			call settst(tst5)
		endif
		if(il==4) then
			stage=2
			level=0.1
			call loadlam
			call settst(tst10)
			level=0.01
			call loadlam
			call settst(tst1)
		endif
	endif

	level=levellist(il)
	levelfudge=levelfudgelist(il)
	
	write(*,*),"XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
	write(*,*),"XXXXXXXXXXXX entering level=", level, "XXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
	write(*,*),"XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
	
	stage=1
	call setcvs
	call setcred
	call preppairs

	if(stage==0) then
		call chksizeopt
		stop
	endif
	
	call setlam0
	call setlam(0)
	print *,"done with stage1"
	

	stage=2
	call setcred
	call preppairs
	call setlam0
	call setlam(0)
	print *,"done with stage2"
	call settst(sstst)
	do in=1,size(nlist)
		n_sample=nlist(in)
		call evalss
	enddo

	if(il<=2)	call mkMCtable

	enddo
	
	call printtime

!	stop if no tables other than baseline table is required for current (k,n0) combo
!	stop			
	
!* now computing tables
	call hmdaloaddata
	
	stage=2
	level=0.05
	call loadlam
	call settst(tst5)

	do il=1,2
		level=levellist(il)
		call loadlam
		stage=2
		call settst(sstst)
		call mkMC2Stable
		call mkMCclustertable
		call saveempMC
		call saveempMCtwosamples
		call saveclustempMC
	enddo

! following routine collects output files and generates Table 2
! all relevant (k,n0) combinations must have been computed already
!	call mkMCktable	
end program